import random

import sympy, numpy
import math
import matplotlib.pyplot as pl
from mpl_toolkits.mplot3d import Axes3D as ax3
import numpy as np
from numpy.linalg import LinAlgError

'''
6. 对于随机生成的(P, q, r, x 0 )，编程实现“梯度下降+精确直线搜索”求解上述优化问题
'''


def getP(dim):
    while 1:
        A = np.random.rand(dim, dim)
        B = np.dot(A, A.transpose())
        C = B + B.T  # makesure symmetric
        try:
            # test whether C is definite
            np.linalg.cholesky(C)  # if there is no error, C is definite
            return C
        except ValueError:
            pass
        except LinAlgError:
            pass


# 最速下降法，二维实验
def SD(x0, G, b, c, N, E):
    f = lambda x: 0.5 * (numpy.dot(numpy.dot(x.T, G), x)) + numpy.dot(b.T, x) + c  # 定义原函数
    f_d = lambda x: numpy.dot(G, x) + b  # 导数
    X = x0
    Y = []  # 存储y的值
    Y_d = []  # 存储y的导数值
    xx = sympy.symarray('xx', (2, 1))
    n = 1
    ee = f_d(x0)

    e = (ee[0] ** 2 + ee[1] ** 2) ** 0.5  # 评价迭代精度
    Y.append(f(x0)[0, 0])
    Y_d.append(e)
    a = sympy.Symbol('a', real=True)  # 最速下降法中的namita
    print('第%d次迭代：e=%d' % (n, e))
    while n < N and e > E:
        n = n + 1
        yy = f(x0 - a * f_d(x0))
        yy_d = sympy.diff(yy[0, 0], a, 1)
        a0 = sympy.solve(yy_d)
        x0 = x0 - a0 * f_d(x0)  # 更新迭代
        X = numpy.c_[X, x0]
        Y.append(f(x0)[0, 0])
        ee = f_d(x0)
        e = math.pow(math.pow(ee[0, 0], 2) + math.pow(ee[1, 0], 2), 0.5)
        Y_d.append(e)
        print(ee)
        print('第%d次迭代：e=%s，lr=%s' % (n, e,str(a0)))
    return X, Y, Y_d


if __name__ == '__main__':
    G = numpy.array([[2, -2], [-2, 4]])
    b = numpy.array([[-4], [0]])
    c = 0
    f = lambda x: 0.5 * (numpy.dot(numpy.dot(x.T, G), x)) + numpy.dot(b.T, x) + c

    f_d = lambda x: numpy.dot(G, x) + b

    x0 = numpy.array([[1], [1]])
    print(f_d(x0))
    N = 40
    E = 10 ** (-6)
    X, Y, Y_d = SD(x0, G, b, c, N, E)
    # 画图
    X = numpy.array(X)
    Y = numpy.array(Y)
    Y_d = numpy.array(Y_d)
    figure1 = pl.figure('trend')
    n = len(Y)
    x = numpy.arange(1, n + 1)
    pl.subplot(2, 1, 1)
    pl.semilogy(x, Y, 'r*')
    pl.xlabel('n')
    pl.ylabel('f(x)')
    pl.subplot(2, 1, 2)
    pl.semilogy(x, Y_d, 'g*')
    pl.xlabel('n')
    pl.ylabel("|f'(x)|")
    figure2 = pl.figure('behave')
    x = numpy.arange(-110, 110, 1)
    y = x
    [xx, yy] = numpy.meshgrid(x, y)
    zz = numpy.zeros(xx.shape)
    n = xx.shape[0]
    for i in range(n):
        for j in range(n):
            xxx = numpy.array([xx[i, j], yy[i, j]])
            zz[i, j] = f(xxx.T)
    ax = ax3(figure2)
    ax.contour3D(xx, yy, zz)
    ax.plot3D(X[0, :], X[1, :], Y, 'ro--')
    pl.show()
